YAPS - Yet Another Positioning Solver

ETN Training sessions 2022, ETN Training Session 2022, České Budějovice, Czechia
Part 3 - Tracking pike in a large lake

Technical University of Denmark (DTU)
hba@aqua.dtu.dk

Setup and needed libraries

Set timezone to UTC and load needed libraries

rm(list=ls())
Sys.setenv(TZ='UTC')
library(data.table)
library(leaflet)
library(sp)
library(lubridate)
library(ggplot2)
library(caTools)
library(viridis)
library(yaps)
library(yapsdata)
source('yaps_ETN_helperFuncs.R')

Description

These data are included in the yaps package (run ?hald for info). Data were collected as part of a larger study including brown trout (smolt and adult), pike and eel in the Danish lake Hald. Lake Hald is large and deep (by Danish standards; 3.42 km²; max depth 34 m), so depth is relevant in order to get best possible track estimation (but ignored here for simplicity). A total of 70 hydrophones (Thelma TBR700) were deployed to cover the entire lake. The data set contains a single week of detections including a single tow track and tracks of a few pike (Esox lucius).

This example focus on:

• Tag 5138: tow track, fixed BI (10 s), gps data available
• Tag 1335: pike, random BI (10-30 s), known sequence
• Tag 1315: pike, random BI (60-120 s), known sequence

An appropriate sync_model is included in the package and will be used for this demonstration to save time. Code to obtain this sync_model can be found in the appendix below.

A look at the data

Notice the object burst_seqs - these are vectors of intervals between consequtive pings from the two transmitters included in this example data. Such data is not available from all manufacturers.
data.tables hydros and detections are as usual.
The object sync_model contains the already optimized synchronization model applicable for these data.

names(hald)

## [1] "hydros"     "detections" "sync_model" "fish"       "burst_seqs"

hydros <- hald$hydros
detections <- hald$detections
sync_model <- hald$sync_model
str(hald$burst_seqs)

## List of 2
##  $ seq1315: int [1:60290] 73 106 110 102 113 89 65 115 103 114 ...
##  $ seq1335: int [1:270460] 14 22 29 18 25 10 24 14 28 15 ...

The study site

Study site geometry is rather complex. .
Note, how a sequential synchronization process in this case could lead to error propagation. Also note, there are lots of impossible sync tag - hydro combinations.

A look at the sync model

There might be room for improvement, but it will suffice for now. Additionally, the hope for this study was not to achieve sub-meter accuracy and precision given the large extent of the study area.

plotSyncModelResids(sync_model, by='overall')

plotSyncModelResids(sync_model, by='quantiles')

plotSyncModelResids(sync_model, by='sync_tag')

plotSyncModelResids(sync_model, by='hydro')

plotSyncModelResids(sync_model, by = 'temporal_hydro')

plotSyncModelResids(sync_model, by = 'temporal_sync_tag')

# # # The following checks might throw an error - looking into it...
plotSyncModelCheck(sync_model, by = "hydro")

plotSyncModelCheck(sync_model, by = "sync_tag")

plotSyncModelCheck(sync_model, by = "sync_bin_sync")

plotSyncModelCheck(sync_model, by = "sync_bin_hydro")

plotSyncModelCheck(sync_model, by = "sync_bin_sync_smooth")

plotSyncModelCheck(sync_model, by = "sync_bin_hydro_smooth")

As part of fitting the synchronization model, we estimated positions of two hydros as their gps-obtained positions were uncertain. We will use these estimated positions for the track estimation.

which(sync_model$inp_synced$dat_tmb_sync$fixed_hydros_vec !=1)

## [1] 32 70

hydros_yaps <- data.table::data.table(sync_model$pl$TRUE_H)
colnames(hydros_yaps) <- c('hx','hy','hz')

cbind(hydros[c(32, 70), c('x','y')], hydros_yaps[c(32,70), c('hx','hy')])

Applying the sync model to all detections

This is easily done, since we already have the sync_model

detections_synced <- applySync(toa=detections, hydros=hydros, sync_model=sync_model)

Estimating tracks using yaps

Tag 5138 - tow track - fixed burst interval

Compiling the data

Remember this is a fixed burst interval transmitter (BI = 10 s).

Let’s have a look at the data. Note, this plot only makes sense for transmitters where the burst interval is set at whole seconds (random or not).

dat_tow <- detections_synced[tag==5138]
gps <- fread("data/hald_gps_test_track.csv")

with(dat_tow, table(hydro_idx))

## hydro_idx
##   1   2   3   4   5   6   7   8   9  53  54  55  56  57  58  59  60  61  62  63 
## 121 241 277 404 420 394 400 242   9  28  73 287 399 402 318  54 206 419 401 202 
##  64  69  70 
##   4 417 409

ggplot(data=dat_tow) + geom_point(aes(x=eposync, y=eposync-floor(eposync), col=factor(hydro_idx)))

Get TOA matrix for YAPS and have a look

toa_tow <- getToaYaps(synced_dat=dat_tow, hydros=hydros, rbi_min=10, rbi_max=10)

## WARNING: pingType not specified in getToaYaps() - will assume 'rbi'. This will become a fatal error in later versions.

plotToa(toa_tow)

plotNobs(toa_tow)

Get inp and run YAPS

set.seed(1)
inp_tow <- getInp(toa=toa_tow, hydros=hydros_yaps, E_dist="Mixture", n_ss=5, pingType="sbi", sdInits=1, ss_data_what="est", ss_data=0, biTable=NULL)
yaps_tow <- runYaps(inp_tow, maxIter=1000, getPlsd=TRUE, getRep=TRUE, silent=TRUE)

## [1] "Pre-flight checkInp() passed!"
## [1] "Running yaps..."
## [1] "...yaps converged (obj: -29506.1936027368) with message: false convergence (8)"

Results

plotYaps(yaps_tow)
lines(Y~X, data=gps, lty=2)

# zoom in 
plotYaps(yaps_tow, xlim=c(520500, 521500), ylim=c(6247400, 6247800))
lines(Y~X, data=gps, lty=2)

plotYaps(yaps_tow, type="coord_X")

plotYaps(yaps_tow, type="coord_Y")

# standard error of position estimates are hardly visible
plot(yaps_tow$plsd$X, type="l")

plot(yaps_tow$plsd$Y, type="l")

# Fixed burst interval - how fixed? Note y-axis is in meter
bi <- diff(yaps_tow$pl$top)
plot((bi - mean(bi)) * 1460, ylab="BI variation (m)", type="h")

# estimated speed of sound - sanity check
yaps_tow$pl$ss

## [1] 1463.449 1462.667 1463.732 1465.495 1464.731

getSS(temp=c(15, 16, 17))

## [1] 1464.915 1467.978 1470.967

# look at the residual
matplot(t((yaps_tow$rep$mu_toa - inp_tow$datTmb$toa))*1460, ylab="TOA residuals (m)")

## Warning in matplot(t((yaps_tow$rep$mu_toa - inp_tow$datTmb$toa)) * 1460, :
## default 'pch' is smaller than number of columns and hence recycled

Tag 5138 tow track - downsampling max_hydro <- 3

Looks good, but this data set is highly saturated! To showcase the benefit of fixed burst interval, we downsample and re-run yaps

plotNobs(toa_tow)

max_hydro <- 3
toa_tow3<- downSampleToa(toa_tow, max_hydro)

set.seed(1)
inp_tow3 <- getInp(toa=toa_tow3, hydros=hydros_yaps, E_dist="Mixture", n_ss=5, pingType="sbi", sdInits=1, ss_data_what="est", ss_data=0, biTable=NULL)
yaps_tow3 <- runYaps(inp_tow3, maxIter=500, getPlsd=TRUE, getRep=TRUE, silent=TRUE)

## [1] "Pre-flight checkInp() passed!"
## [1] "Running yaps..."
## [1] "...yaps converged (obj: -5706.29991463201) with message: relative convergence (4)"

plotYaps(yaps_tow3, xlim=c(520500, 521500), ylim=c(6247400, 6247800))
lines(Y~X, data=gps, lty=2)

plotYaps(yaps_tow3, type="coord_X")

plotYaps(yaps_tow3, type="coord_Y")

Tag 5138 tow track - downsampling max_hydro <- 1

For fixed burst interval transmitters, we can go really low on detections.

plotNobs(toa_tow)

max_hydro <- 1
toa_tow1<- downSampleToa(toa_tow, max_hydro)

set.seed(1)
inp_tow1 <- getInp(toa=toa_tow1, hydros=hydros_yaps, E_dist="Mixture", n_ss=5, pingType="sbi", sdInits=1, ss_data_what="est", ss_data=0, biTable=NULL)
yaps_tow1 <- runYaps(inp_tow1, maxIter=500, getPlsd=TRUE, getRep=TRUE, silent=TRUE)

## [1] "Pre-flight checkInp() passed!"
## [1] "Running yaps..."
## [1] "...yaps converged (obj: -1669.02009069074) with message: false convergence (8)"

plotYaps(yaps_tow1, xlim=c(520500, 521500), ylim=c(6247400, 6247800))
lines(Y~X, data=gps, lty=2)

plotYaps(yaps_tow1, type="coord_X")

plotYaps(yaps_tow1, type="coord_Y")

Comparing standard errors of position estimates

par(mfrow=c(2,2))
plotYaps(yaps_tow, type="coord_X", main="X all data")
plotYaps(yaps_tow, type="coord_Y", main="Y all data")
plotYaps(yaps_tow1, type="coord_X", main="X max_hydro=1")
plotYaps(yaps_tow1, type="coord_Y", main="Y max_hydro=1")

par(mfrow=c(1,1))

par(mfrow=c(3,1))
ylim <- c(0, 5)
plot(yaps_tow$plsd$X,  type="h", ylim=ylim, main="Full")
plot(yaps_tow3$plsd$X, type="h", ylim=ylim, main="max_hydro=3")
plot(yaps_tow1$plsd$X, type="h", ylim=ylim, main="max_hydro=1")

par(mfrow=c(1,1))

Note on the downsampling

One reason we can estimate a decent track even with max_hydro = 1 is, that multiple hydrophones are detecting the signals. It will, of course, not work if there only was a single hydrophone in the water.

Tag 1335 - pike - random BI sequence - assuming unknown sequence

Remember this transmitter is a random burst interval (BI 10 - 30 s), but we know the burst sequence. However, first we assume, we don’t and treat it like a random BI transmitter with unknown burst sequence.

Compiling the data

dat1335 <- detections_synced[tag==1335]
ping_type <- 'rbi'
rbi_min <- 10
rbi_max <- 30

# get number of detections per hydro per hour
dat1335[, hour := floor_date(ts, unit="hour")]
summ <- dat1335[, .N, by=c('hydro_idx', 'hour')]
ggplot(data=summ) + geom_tile(aes(x=hour, y=hydro_idx, fill=N)) + scale_fill_viridis()

# Get TOA matrix for YAPS
toa1335_rbi <- getToaYaps(synced_dat=dat1335, hydros=hydros_yaps, rbi_min=rbi_min, rbi_max=rbi_max)

## WARNING: pingType not specified in getToaYaps() - will assume 'rbi'. This will become a fatal error in later versions.

# Have a look at the raw data...
plotToa(toa1335_rbi)

plotNobs(toa1335_rbi)

More than 30.000 pings will take pretty long time to process, so we use a small a subset of the data.

ping_1 <- 501
n_ping <- 1000
toa1335_rbi_yaps <- toa1335_rbi[ping_1:(ping_1 + n_ping - 1), ]
plotToa(toa1335_rbi_yaps)

plotNobs(toa1335_rbi_yaps)

Get inp and run YAPS

set.seed(1)
ping_type <- 'rbi'
inp1335_rbi <- getInp(toa=toa1335_rbi_yaps, hydros=hydros_yaps, E_dist="t", n_ss=5, pingType=ping_type, rbi_min=rbi_min, rbi_max=rbi_max, sdInits=1, ss_data_what="est", ss_data=0, biTable=NULL)

# We cheat and use fixed initial values to ensure convergence for this specific example. 
# DO NOT DO THIS when analysing own data!!!
inp1335_rbi$inits <- c(-0.5842236,  0.5724419, -0.8182362,  4.9533368)
yaps1335_rbi <- runYaps(inp1335_rbi, maxIter=1000, getPlsd=TRUE, getRep=TRUE, silent=TRUE)

## [1] "Pre-flight checkInp() passed!"
## [1] "Running yaps..."
## [1] "...yaps converged (obj: -11416.6759567938) with message: false convergence (8)"

Results

# basic plotting
plotYaps(yaps1335_rbi)

# zoom in 
plotYaps(yaps1335_rbi, xlim=c(520500, 521500), ylim=c(6247400, 6247800))

plotYaps(yaps1335_rbi, type="coord_X")

plotYaps(yaps1335_rbi, type="coord_Y")

# uncertainties... standard error of position estimates
plot(yaps1335_rbi$plsd$X, type="h")

plot(yaps1335_rbi$plsd$Y, type="h")

# estimated speed of sound - sanity check
yaps1335_rbi$pl$ss

## [1] 1470.552 1470.477 1460.903 1461.303 1460.754

getSS(temp=c(15, 17))

## [1] 1464.915 1470.967

# look at the residual
matplot(t((yaps1335_rbi$rep$mu_toa - inp1335_rbi$datTmb$toa))*1460, ylab="TOA residuals (m)")

## Warning in matplot(t((yaps1335_rbi$rep$mu_toa - inp1335_rbi$datTmb$toa)) * :
## default 'pch' is smaller than number of columns and hence recycled

Tag 1335 - pike - random BI - utilizing known sequence

To utilize the known sequence of burst intervals, we need to align the observed data with the sequence. For this, we can use the function alignBurstSeq() from yaps. We can also take advantage of knowing the burst sequence when compiling the TOA-marix and we need extra info to feed into getInp().

Compiling the data

dat1335 <- detections_synced[tag==1335]
ping_type <- 'pbi'
rbi_min <- 10
rbi_max <- 30

seq1335 <- hald$burst_seqs$seq1335 # just a long sequence of random numbers...
length(seq1335)

## [1] 270460

dat1335_kbi <- alignBurstSeq(synced_dat=dat1335, burst_seq=seq1335, seq_lng_min=25, rbi_min=rbi_min, rbi_max=rbi_max, plot_diag=TRUE)

toa1335_kbi_list <- buildToaKnownSeq(seq=seq1335, aligned_dat=dat1335_kbi, hydros=hydros_yaps)

toa1335_kbi <- toa1335_kbi_list$toa
seq1335_kbi <- toa1335_kbi_list$seq

More than 30 K pings will still take a long time to finish - we use the same subset as above.

ping_1 <- 501
n_ping <- 1000
toa1335_kbi_yaps <- toa1335_kbi[ping_1:(ping_1 + n_ping - 1), ]
biTable1335 <- seq1335_kbi[ping_1:(ping_1 + n_ping - 1) ]
plotToa(toa1335_kbi_yaps)

plotNobs(toa1335_kbi_yaps)

Get inp and run YAPS

Note, that ping_type now is 'pbi'and we include the BItable in getInp(). Everything else is the same as above.

set.seed(1)
inp1335_kbi <- getInp(toa=toa1335_kbi_yaps, hydros=hydros_yaps, E_dist="t", n_ss=5, pingType=ping_type, rbi_min=rbi_min, rbi_max=rbi_max, sdInits=1, ss_data_what="est", ss_data=0, biTable=biTable1335)

yaps1335_kbi <- runYaps(inp1335_kbi, maxIter=1000, getPlsd=TRUE, getRep=TRUE, silent=TRUE)

## [1] "Pre-flight checkInp() passed!"
## [1] "Running yaps..."
## [1] "...yaps converged (obj: -20284.7251811341) with message: relative convergence (4)"

Results

Quick look at the results

# basic plotting
plotYaps(yaps1335_kbi)

# zoom in 
plotYaps(yaps1335_kbi, xlim=c(520500, 521500), ylim=c(6247400, 6247800))

plotYaps(yaps1335_kbi, type="coord_X")

plotYaps(yaps1335_kbi, type="coord_Y")

# uncertainties... standard error of position estimates
plot(yaps1335_kbi$plsd$X, type="h")

plot(yaps1335_kbi$plsd$Y, type="h")

# estimated speed of sound - sanity check
yaps1335_kbi$pl$ss

## [1] 1469.958 1470.979 1458.072 1461.020 1460.114

getSS(temp=c(15, 17))

## [1] 1464.915 1470.967

# look at the residual
matplot(t((yaps1335_kbi$rep$mu_toa - inp1335_kbi$datTmb$toa))*1460, ylab="TOA residuals (m)")

## Warning in matplot(t((yaps1335_kbi$rep$mu_toa - inp1335_kbi$datTmb$toa)) * :
## default 'pch' is smaller than number of columns and hence recycled

Tag 1335 pike - downsampling max_hydro <- 3

Knowing the burst interval sequence, gives opportunity to estimate tracks even when data density is low.

set.seed(1)
ping_type <- 'pbi'

plotNobs(toa1335_kbi_yaps)

max_hydro <- 3
toa1335_kbi_down3 <- downSampleToa(toa1335_kbi_yaps, max_hydro)

inp1335_kbi_down3 <- getInp(toa=toa1335_kbi_down3, hydros=hydros_yaps, E_dist="t", n_ss=5, pingType=ping_type, rbi_min=rbi_min, rbi_max=rbi_max, sdInits=1, ss_data_what="est", ss_data=0, biTable=biTable1335)

yaps1335_kbi_down3 <- runYaps(inp1335_kbi_down3, maxIter=500, getPlsd=TRUE, getRep=TRUE, silent=TRUE)

## [1] "Pre-flight checkInp() passed!"
## [1] "Running yaps..."
## [1] "...yaps converged (obj: -11913.2349725942) with message: relative convergence (4)"

plotYaps(yaps1335_kbi_down3)

# zoom in 
plotYaps(yaps1335_kbi_down3, xlim=c(520500, 521500), ylim=c(6247400, 6247800))

plotYaps(yaps1335_kbi_down3, type="coord_X")

plotYaps(yaps1335_kbi_down3, type="coord_Y")

Tag 1335 pike - downsampling max_hydro <- 1

Trying at very low data density

set.seed(1)
ping_type <- 'pbi'

plotNobs(toa1335_kbi_yaps)

max_hydro <- 1
toa1335_kbi_down1 <- downSampleToa(toa1335_kbi_yaps, max_hydro)

inp1335_kbi_down1 <- getInp(toa=toa1335_kbi_down1, hydros=hydros_yaps, E_dist="Mixture", n_ss=5, pingType=ping_type, rbi_min=rbi_min, rbi_max=rbi_max, sdInits=1, ss_data_what="est", ss_data=0, biTable=biTable1335)

yaps1335_kbi_down1 <- runYaps(inp1335_kbi_down1, maxIter=500, getPlsd=TRUE, getRep=TRUE, silent=TRUE)

## [1] "Pre-flight checkInp() passed!"
## [1] "Running yaps..."
## [1] "...yaps converged (obj: -3645.53522305893) with message: false convergence (8)"

# zoom in 
plotYaps(yaps1335_kbi_down1, xlim=c(520500, 521500), ylim=c(6247400, 6247800))

plotYaps(yaps1335_kbi_down1, type="coord_X")

plotYaps(yaps1335_kbi_down1, type="coord_Y")

Comparing unknown vs known BI sequences

par(mfrow=c(2,3))
ylim=c(0,5)
plot(yaps1335_kbi$plsd$X, type="h", ylim=ylim, main="SD all")
plot(yaps1335_kbi_down1$plsd$X, type="h", ylim=ylim, main = "SD nObs <= 1")
plot(yaps1335_rbi$plsd$X, type="h", ylim=ylim, main="SD rbi")
plot(yaps1335_kbi$plsd$Y, type="h", ylim=ylim)
plot(yaps1335_kbi_down1$plsd$Y, type="h", ylim=ylim)
plot(yaps1335_rbi$plsd$Y, type="h", ylim=ylim)

par(mfrow=c(1,1))

Tag 1335 - the full track

The full track can be estimated using this code:

set.seed(42)
ping_type <- 'pbi'
rbi_min <- 10
rbi_max <- 30
seq1335 <- hald$burst_seqs$seq1335 # just a long sequence of random numbers...
dat1335_kbi <- alignBurstSeq(synced_dat=dat1335, burst_seq=seq1335, seq_lng_min=10, rbi_min=rbi_min, rbi_max=rbi_max, plot_diag=TRUE)
toa1335_kbi_list <- buildToaKnownSeq(seq=seq1335, aligned_dat=dat1335_kbi, hydros=hydros_yaps)
toa1335_kbi <- toa1335_kbi_list$toa
seq1335_kbi <- toa1335_kbi_list$seq
biTable1335 <- seq1335_kbi

ping_1 <- 1
n_ping <- nrow(toa1335_kbi)
toa1335_kbi_yaps <- toa1335_kbi[ping_1:(ping_1 + n_ping - 1), ]
biTable1335 <- seq1335_kbi[ping_1:(ping_1 + n_ping - 1) ]
plotToa(toa1335_kbi_yaps)
plotNobs(toa1335_kbi_yaps)

inp1335_kbi <- getInp(toa=toa1335_kbi_yaps, hydros=hydros_yaps, E_dist="t", n_ss=5, pingType=ping_type, rbi_min=rbi_min, rbi_max=rbi_max, sdInits=1, ss_data_what="est", ss_data=0, biTable=biTable1335)
#inp1335_kbi$inits[1] <- 1
yaps1335_kbi <- runYaps(inp1335_kbi, maxIter=1000, getPlsd=TRUE, getRep=TRUE, silent=FALSE)

Lazy people can just grab it here:

load('data/yaps1335_kbi_fullTrack.rObj')
plotYaps(yaps1335_kbi)

plotYaps(yaps1335_kbi, type="coord_X")

plotYaps(yaps1335_kbi, type="coord_Y")

Tag 1315 - pike - random BI - assuming unknown sequence

This transmitter is a random BI 60 - 120 s.

dat1315 <- detections_synced[tag==1315]
rbi_min <- 60
rbi_max <- 120
dat1315[, hour := floor_date(ts, unit="hour")]
summ <- dat1315[, .N, by=c('hydro_idx', 'hour')]
ggplot(data=summ) + geom_tile(aes(x=hour, y=hydro_idx, fill=N)) + scale_fill_viridis()

# First, assume burst sequence is unknown
# Get TOA matrix for YAPS
ping_type <- 'rbi'
toa1315_rbi <- getToaYaps(synced_dat=dat1315, hydros=hydros_yaps, rbi_min=rbi_min, rbi_max=rbi_max)

## WARNING: pingType not specified in getToaYaps() - will assume 'rbi'. This will become a fatal error in later versions.

# Have a look at the raw data...
plotToa(toa1315_rbi)

plotNobs(toa1315_rbi)

nobs <- getNobs(toa1315_rbi)
sum(nobs>=3) / length(nobs)

## [1] 0.6995993

# < 10.000 pings is doable, but we will sub sample to save time
set.seed(1)
ping_1 <- 1
n_ping <- 1000
toa1315_rbi_yaps <- toa1315_rbi[ping_1:(ping_1 + n_ping - 1), ]
plotToa(toa1315_rbi_yaps)

plotNobs(toa1315_rbi_yaps)

nobs <- getNobs(toa1315_rbi_yaps)
sum(nobs>=3) / length(nobs)

## [1] 0.46

# a bit more challenging than 1335

# get inp to feed into YAPS and run YAPS
inp1315_rbi <- getInp(toa=toa1315_rbi_yaps, hydros=hydros_yaps, E_dist="t", n_ss=5, pingType=ping_type, rbi_min=rbi_min, rbi_max=rbi_max, sdInits=1, ss_data_what="est", ss_data=0, biTable=NULL)
# inp1315_rbi$inits[1] <- 1
yaps1315_rbi <- runYaps(inp1315_rbi, maxIter=1000, getPlsd=TRUE, getRep=TRUE, silent=TRUE)

## [1] "WARNING: inp$datTmb$rbi_max < max(diff(inp$params$top)) | 126 < 130.123713254929"
## [1] "Pre-flight checkInp() passed!"
## [1] "Running yaps..."
## [1] "...yaps converged (obj: -1359.96122380651) with message: false convergence (8)"

# basic plotting
plotYaps(yaps1315_rbi)

plotYaps(yaps1315_rbi, type="coord_X")

plotYaps(yaps1315_rbi, type="coord_Y")

# uncertainties... standard error of position estimates
plot(yaps1315_rbi$plsd$X, type="l")

plot(yaps1315_rbi$plsd$Y, type="l")

plotNobs(toa1315_rbi_yaps)

# estimated speed of sound - sanity check
yaps1315_rbi$pl$ss

## [1] 1473.040 1473.002 1473.244 1473.812 1473.679

# look at the residual
matplot(t((yaps1315_rbi$rep$mu_toa - inp1315_rbi$datTmb$toa))*1450)

## Warning in matplot(t((yaps1315_rbi$rep$mu_toa - inp1315_rbi$datTmb$toa)) * :
## default 'pch' is smaller than number of columns and hence recycled

Tag 1315 - pike - random BI - utilizing known sequence

ping_type <- 'pbi'
seq1315 <- hald$burst_seqs$seq1315 # just a long sequence of random numbers...

# we need to figure out where our data is on this sequence
dat1315_kbi <- alignBurstSeq(synced_dat=dat1315, burst_seq=seq1315, seq_lng_min=10, rbi_min=rbi_min, rbi_max=rbi_max, plot_diag=TRUE)

toa1315_kbi_list <- buildToaKnownSeq(seq=seq1315, aligned_dat=dat1315_kbi, hydros=hydros_yaps)
str(toa1315_kbi_list)

## List of 2
##  $ toa: num [1:7238, 1:70] NA NA NA NA NA NA NA NA NA NA ...
##  $ seq: int [1:7238] 96 120 87 80 117 118 99 87 78 104 ...

toa1315_kbi <- toa1315_kbi_list$toa
seq1315_kbi <- toa1315_kbi_list$seq
biTable1315 <- seq1315_kbi

# subset data...
ping_1 <- 1
n_ping <- 1000
toa1315_kbi_yaps <- toa1315_kbi[ping_1:(ping_1 + n_ping - 1), ]
biTable1315 <- seq1315_kbi[ping_1:(ping_1 + n_ping - 1) ]
plotToa(toa1315_kbi_yaps)

plotNobs(toa1315_kbi_yaps)

sum(getNobs(toa1315_kbi_yaps) >= 3) / nrow(toa1315_kbi_yaps)

## [1] 0.459

# get inp to feed into YAPS and run YAPS
set.seed(1)
inp1315_kbi <- getInp(toa=toa1315_kbi_yaps, hydros=hydros_yaps, E_dist="t", n_ss=5, pingType=ping_type, rbi_min=rbi_min, rbi_max=rbi_max, sdInits=1, ss_data_what="est", ss_data=0, biTable=biTable1315)
# inp1315_kbi$inits[1] <- 1
yaps1315_kbi <- runYaps(inp1315_kbi, maxIter=1000, getPlsd=TRUE, getRep=TRUE, silent=TRUE)

## [1] "Pre-flight checkInp() passed!"
## [1] "Running yaps..."
## [1] "...yaps converged (obj: -10243.9825865498) with message: relative convergence (4)"

# basic plotting
plotYaps(yaps1315_kbi, type="coord_X")

plotYaps(yaps1315_kbi, type="coord_Y")

plotYaps(yaps1315_kbi)

# uncertainties... standard error of position estimates - compare to ping_type = 'rbi'
par(mfrow=c(2,2))
ylim <- c(0,50)
plot(yaps1315_kbi$plsd$X, type="h", ylim=ylim)
plot(yaps1315_rbi$plsd$X, type="h", ylim=ylim)
plot(yaps1315_kbi$plsd$Y, type="h", ylim=ylim)
plot(yaps1315_rbi$plsd$Y, type="h", ylim=ylim)

par(mfrow=c(1,1))

The full track

The full track can be estimated using this code - will probably take some 5 - 10 minutes.

set.seed(1)
ping_type <- 'pbi'
rbi_min <- 60
rbi_max <- 120
seq1315 <- hald$burst_seqs$seq1315 # just a long sequence of random numbers...
dat1315_kbi <- alignBurstSeq(synced_dat=dat1315, burst_seq=seq1315, seq_lng_min=10, rbi_min=rbi_min, rbi_max=rbi_max, plot_diag=TRUE)
toa1315_kbi_list <- buildToaKnownSeq(seq=seq1315, aligned_dat=dat1315_kbi, hydros=hydros_yaps)
toa1315_kbi <- toa1315_kbi_list$toa
seq1315_kbi <- toa1315_kbi_list$seq
biTable1315 <- seq1315_kbi

ping_1 <- 1
n_ping <- nrow(toa1315_kbi)
toa1315_kbi_yaps <- toa1315_kbi[ping_1:(ping_1 + n_ping - 1), ]
biTable1315 <- seq1315_kbi[ping_1:(ping_1 + n_ping - 1) ]
plotToa(toa1315_kbi_yaps)
plotNobs(toa1315_kbi_yaps)

inp1315_kbi <- getInp(toa=toa1315_kbi_yaps, hydros=hydros_yaps, E_dist="Mixture", n_ss=5, pingType=ping_type, rbi_min=rbi_min, rbi_max=rbi_max, sdInits=1, ss_data_what="est", ss_data=0, biTable=biTable1315)
# inp1315_kbi$inits[1] <- 1
yaps1315_kbi <- runYaps(inp1315_kbi, maxIter=1000, getPlsd=TRUE, getRep=TRUE, silent=TRUE)

Lazy people can just grab it here:

load('data/yaps1315_kbi_fullTrack.rObj')
plotYaps(yaps1315_kbi)

plotYaps(yaps1315_kbi, type="coord_X")

plotYaps(yaps1315_kbi, type="coord_Y")

Appendix

Obtaining the sync_model for the hald data

If you try to sync these data (which is encouraged), be aware of impossible sync tag - hydro combos. Also note that positions of hydro idx 32 and 70 are uncertain, so they should be estimated, i.e. removed from the fixed_hydros_idx vector. These settings should give a good sync model

sync_dat <- list(detections = hald$detections, hydros = hald$hydros)
    

max_epo_diff <- 30
min_hydros <- 5
time_keeper_idx <- 1
fixed_hydros_idx <- 1:70
fixed_hydros_idx <- fixed_hydros_idx[!fixed_hydros_idx %in% c(32, 70)]
n_offset_day <- 2
n_ss_day <- 1
inp_sync <- getInpSync(sync_dat, max_epo_diff, min_hydros, time_keeper_idx, fixed_hydros_idx, n_offset_day, n_ss_day, keep_rate=25)
sync_model0 <- getSyncModel(inp_sync, silent=TRUE)
sync_model1 <- fineTuneSyncModel(sync_model0, eps_threshold = 100)
    
# Plot model residuals and model check plots to ensure the synchronization process was successful...
plotSyncModelResids(sync_model, by = 'overall')    
plotSyncModelResids(sync_model, by = 'quantiles')
plotSyncModelResids(sync_model, by = 'sync_tag')      
plotSyncModelResids(sync_model, by = 'hydro')  
plotSyncModelResids(sync_model, by = 'temporal_hydro')
plotSyncModelResids(sync_model, by = 'temporal_sync_tag')

knitr::purl(‘part3_hald.Rmd’, documentation=2)